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ABSTRACT 

Dynamical instability is shown to occur in differentially rotating polytropes with N — 3.33 
and T/|W| > 0.14. This instability has a strong m = 1 mode, although the m — 2,3, and 4 
modes also appear. Such instability may allow a centrifugally-hung core to begin collapsing to 
neutron star densities on a dynamical timescale. The gravitational radiation emitted by such 
unstable cores may be detectable with advanced ground-based detectors, such as LIGO II. If the 
instability occurs in a supermassive star, it may produce gravitational radiation detectable by 
the space-based detector LISA. 

Subject headings: instabilities — gravitation — hydrodynamics — stars: rotation — stars: neutron 



1. Introduction 

Rotational instabilities are potentially impor- 
tant in the evolution of massive stellar cores. For 
example, the core of a massive star that has been 
prevented from collapsing to neutron star densities 
by centrifugal forces may rotate rapidly enough 
for rotational instabilities to develop (Hayashi, 
Eriguchi, & Hashimoto 1998, 1999). Massive 
cores spun up by accretion from a binary compan- 
ion (Wagoner 1984) and the remnants of compact 
binary mergers (Rasio & Shapiro 1992; Zhuge, 
Centrella, & McMillan 1994) may also reach fast 
enough rotation rates to become unstable. Rota- 
tional instabilities may produce detectable gravi- 
tational radiation (Schutz 1989; Thorne 1996). 
If enough angular momentum is shed, full collapse 
to a neutron star or black hole may occur. 

We focus here on global rotational instabilities 
that arise in fluids from growing azimuthal modes 
e im<t> (Tassoul 1978). Dynamical instabilities are 
driven by hydrodynamics and gravity, and grow 
on the order of the dynamical timescales of the 
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system. In objects with Maclaurin-like rotation 
laws, dynamical rotational instabilities are gen- 
erally expected to arise at fairly high values of 
(3 = T/\W\ > 0.27 in Newtonian gravity (see 
New, Centrella, & Tohline and references therein), 
and f3 > 0.25 — 0.26 in general relativity (Shibata, 
Baumgarte, & Shapiro 2000). Here, T is the ro- 
tational kinetic energy, and W is the gravitational 
potential energy. Dynamical instability may set 
in at lower values of j3 for thick, self-gravitating 
disks (Pickett, Durisen, & Davis 1996; Woodward, 
Tohline, & Hachisu 1994) and in tori (Tohline & 
Hachisu 1990). Secular instabilities develop on 
longer timescales, and arise from dissipative pro- 
cesses such as viscosity and gravitational radiation 
reaction. They are expected to set in for relatively 
low values of j3 > 0.14, although the rotation law, 
polytropic index, and general relativity can affect 
the value of the instability limit (Imamura, et al. 
1995; Stergioulas & Friedman 1998). 

Studies of scenarios for the formation of centrifugally- 
hung cores have generally concluded that these ob- 
jects will have f3 < 0.27 (Tohline 1984; Eriguchi 
& Miiller 1985; Midler & Eriguchi 1985). Direct 
numerical simulations of axisymmetric stellar col- 
lapse by Zwerger & Miiller (1997) also indicate 
that, for the cases studied, it is difficult for collaps- 
ing cores to reach and exceed f3 ~ 0.27. Overall, it 
is generally assumed that centrifugally-hung cores 
will evolve secularly, as in the "fizzler" scenario 
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(Hayashi, Eriguchi, & Hashimoto 1998, 1999). 

In this Letter, we report on new 3-D numeri- 
cal simulations of rotating stellar cores that ex- 
hibit dynamical rotational instability at relatively 
low values of f3 > 0.14. These models are com- 
puted using Newtonian gravity with no backre- 
action; this is a reasonable approximation for a 
stellar core of mass M ~ 1.4M Q hung up at a ra- 
dius R > 100km, so that (GM/Rc 2 ) < 0.02. Our 
results show that compact objects can become dy- 
namically unstable in a region of parameter space 
previously thought to be solely the province of sec- 
ular instabilities. 

2. Numerical Simulations 

Rotating stellar cores are initially assumed 
to be axisymmetric equilibrium polytropes with 
equation of state P = Kp v = Kp 1+1 / N , with 
r = 1.3 (N = 3.33). The angular velocity dis- 
tribution is given by the so-called j— constant ro- 
tation law, Q 2 = j 2 /(d 2 + vj 2 ) 2 , where vj is the 
cylindrical radius and d is an arbitrary constant 
(Hachisu 1986). As d — > 0, the specific angu- 
lar momentum approaches the constant value j . 
Here, we use d — 0.2. The models are computed 
on a uniformly spaced (w, z) grid using the self- 
consistent field method of Hachisu (1986). Vary- 
ing the axis ratio yields equilibria with different 
values of /?; see New, Centrella, & Tohline (2000) 
for details. Figure 1 shows density contours in the 
x — z plane for the 4 models that we evolved here, 
13 = 0.090,0.12,0.14, and 0.18. Notice that, with 
the exception of the f3 = 0.090 model, the density 
maxima are toroidal; similar structures were used 
as initial data in the collapse models of Rampp, 
Ruffert, & Miiller (1998). 

These initial models were introduced into 2 dif- 
ferent 3-D hydrocodes. The first of these is the 
L code discussed in New, Centrella, & Tohline 
(2000), with a cylindrical (w, z, <f>) grid of 64 x 64 x 
128 zones. The second is the PPM code discussed 
by Brown (2000), with a Cartesian (x, y, z) grid of 
128 x 128 x 128 zones. Both codes have the same 
resolution along the x and z— axes, with the equa- 
torial radius Re of the initial model extending out 
to zone 48 in x-direction. The cylindrical code 
imposes equatorial plane symmetry, whereas the 
Cartesian code imposes no symmetries. Random 
perturbations of 1% were imposed on the densi- 



ties, and the models were evolved forward in time. 

3. Results 

The model with initial j3 = 0.14 exhibits a dy- 
namical instability that is very similar in simu- 
lations performed with both the cylindrical and 
Cartesian hydrocodes. These runs were stopped 
when the loss of mass from the edge of the grid be- 
came significant. Figure 2 shows density contours 
in the equatorial plane from the run with the cylin- 
drical code. The green (lightest greyscale) con- 
tour in the first frame delineates the toroidal den- 
sity maximum seen in Fig. 1. As time proceeds, 
this torus "pinches off" to leave a single high den- 
sity region, indicative of an m = 1 mode. At late 
times, the appearance of higher density contours 
signals that the dense region is starting to collapse. 
Time is measured in units of tr> = {Re/GM) 1 / 2 , 
which is the dynamical time for a sphere of mass M 
and radius Re- For Re ~ 200 km, to ~ 6 ms. 3-D 
animations of the evolution of this model, using 
the same contour levels as in Figure 2 and look- 
ing down the z— axis towards the equatorial plane, 
are available for both the cylindrical and Carte- 
sian runs. Additional animations, in which the 
model is viewed from the side, are also available 
for the cylindrical and Cartesian runs. [Note to 
editors: Here, the words "cylindrical" and "Carte- 
sian" should link to the animations in each case.] 

We can quantify the instability by examining 
various Fourier components in the density distri- 
bution. We examine the density in a ring of fixed 
vj and z using a complex azimuthal Fourier de- 
composition, where the amplitudes C m of the var- 
ious components m are defined by C m (vj, z) = 
(1/2tt) f* w p(w,z,ip)e- im vdtp (Tohline, Duriscn, 
& McCollough 1985). The normalized ampli- 
tude is then A m = C m /Co, where Cq{w,z) = 
p(w, z) is the mean density in the ring, and the 
phase angle of the m th component is <f> m (vj, z) = 
tan -1 [lm(A m )/Re(A m )]. We write the phase an- 
gle <j) m — <7 m t, where er m = dcj>/dt is the eigenfrc- 
quency and W rn — cr rn /m is the pattern speed of 
the m th mode. 

Figures 3 and 4 show the growth of the ampli- 
tudes \A m \ for the first four Fourier components in 
the runs with = 0.14 on the cylindrical and 

Cartesian hydrocodes, respectively. For the cylin- 
drical code, these amplitudes were calculated in 
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the equatorial plane in a ring of width Azu = 1 /48 
at radius w = 0.32; see New, Centrella, & Tohline 
(2000) for details. For the Cartesian code, the 
amplitudes were computed on a circle of radius 
w = 0.32 using a nonuniform discretization, which 
avoids grid boundaries, and a linear interpolation 
of density. 

Similar results were found in rings at other val- 
ues of zu. The evolution in both runs is dominated 
by an exponentially growing m — 1 mode (thick 
solid line). Once this mode is well into its expo- 
nential growth phase, an to = 2 mode also devel- 
ops, followed at later times by to = 3 and m = 4 
modes. (The constant amplitude m — 4 signal in 
the Cartesian run is due to the geometry of the nu- 
merical grid.) Since these modes grow rapidly on 
dynamical timescales, they signal dynamical in- 
stability. Numerical values for the growth rates, 
eigenfrequencies, and pattern speeds are shown in 
Table 1 for to = 1 and m = 2. Notice that the 
cylindrical and Cartesian runs yield similar results 
and that the pattern speeds W\ ~ W2, indicating 
that the to = 2 mode is a harmonic of the to = 1 
mode. Longer runs with larger grids and higher 
resolution are needed to obtain reliable values for 
the to = 3 and to = 4 modes. 

The model with = 0.18 also exhibits dynam- 
ical instability. For the evolutions run with both 
codes, the m = 1 mode grew at very close to the 
same rate, starting at about the same time; see Ta- 
ble 1. In the run on the cylindrical code, the modes 
to = 1, 2, 3, and 4 appear sequentially in time and 
W\ ~ W2, as in the model with (3 = 0.14. In the 
run on the Cartesian code, the to = 1 and m = 2 
modes both grew at about the same rate, starting 
at about the same time; however, W\ 7^ W2- At 
this stage, we do not understand this difference in 
the behavior of the to = 2 mode in this case. We 
plan to investigate it using higher resolution sim- 
ulations. The (3 — 0.12 model was run for <~ 40t£,, 
at which time the to = 1 mode was just begin- 
ning to grow; longer runs with higher resolution 
are needed to confirm the development of insta- 
bility in this case. The model with (3 — 0.09 was 
run for ~ SSi^ 1 , and showed no growing modes. 
In the event that we can confirm instability for 
(3 = 0.12, we will return to the (3 = 0.09 case and 
run it further. 

Previous studies of the to = 2 bar-mode insta- 
bility using the cylindrical hydrocode showed that 



significant motion of the system center of mass 
could develop at late times, resulting in a spurious 
to = 1 signal (New, Centrella, & Tohline 2000). 
For the simulations reported here, we monitored 
the position of the overall center of mass and ver- 
ified that it underwent no systematic motion dur- 
ing the development of the to = 1 mode, in both 
the cylindrical and Cartesian codes. When the 
runs with T/\W\ = 0.14 and T/\W\ = 0.18 were 
stopped, some small, spurious center of mass mo- 
tion (< one zone) had developed in both codes. 
This effect is possibly related to the loss of mass 
from the grid as the outer regions expand; runs 
with larger grids are needed to determine the late 
time behavior of these models. Previous stud- 
ies of the bar instability using the Cartesian code 
(Brown 2000) showed that angular momentum 
losses due to numerical inaccuracies can be signf- 
icant. For the (3 — 0.14 and (3 = 0.18 simulations 
with the Cartesian code, the artificial angular mo- 
mentum losses amounted to less than 0.4% and 
1.5%, respectively. 

4. Discussion 

Our results demonstrate that dynamical insta- 
bility can occur in differentially rotating poly- 
tropes with low values of (3 > 0.14, and that this 
instability has a strong, to = 1 character. Figure 1 
shows that the unstable axisymmetric equilibria 
each have a torus of dense material centered on 
the rotation axis within the model. This suggests 
that this instability may be related to the ones 
found by Woodward, Tohline, & Hachisu (1994) 
and Tohline & Hachisu (1990) in their studies 
of polytropic tori with TV = 3/2. We note that 
Pickett, Durisen, & Davis (1996) also found a 
dominant to = 1 instability that set in at a rela- 
tively low (3 > 0.20, in a centrally condensed model 
(without an off-center density maximum). Their 
model was an N = 3/2 polytrope with an n' = 2 
rotation law. The n' = 2 rotation law produces 
configurations with a strong concentration of an- 
gular momentum in the outer regions of the model. 
We plan to carry out more detailed studies to in- 
vestigate the character of the unstable modes seen 
in our simulations and their properties for vari- 
ous values of TV and the j-constant rotation law 
parameter d. 

If such instability occurs in a centrifugally-hung 
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core, collapse to neutron star densities may result. 
Our simulations do show that the density is in- 
creasing at the end of the unstable runs. Further 
studies are needed to see how dense the remnant 
actually becomes, and whether the m = 1 mode 
will result in the dense region moving at a veloc- 
ity comparable to those of actual neutron stars 
(Popov, S. ct al. 2000). 

Dynamical instability will also produce gravita- 
tional radiation. Although longer runs are needed 
to obtain the full gravitational waveforms, we can 
estimate the properties of the gravitational wave 
emission from the initial stages of the instability. 
For a massive stellar core with M ~ 1.4M and 
Re ~ 200km, the peak amplitude will occur at a 
frequency of roughly / ~ 200Hz. This peak am- 
plitude will be h ~ lO" 24 ^ 1 for = 0.14, and 
h ~ lO -23 ^ 1 for = 0.18. Here, r 20 is the dis- 
tance to the source in units of 20 Mpc. Emis- 
sion from such unstable cores may be detectable 
with advanced ground-based interferometers, such 
as LIGO-II. An even more optimistic scenario for 
detectable gravitational radiation occurs if these 
instabilities arise in supermassive stars 5 (New & 
Shapiro 2000). For example, if the instability oc- 
curs when a supermassive star contracts to the 
point that (GM/Rc 2 ) ~ 1/15, an approximate 
value for uniformly rotating stars (Baumgarte & 
Shapiro 1999), we estimate the frequency of the 
gravitational radiation to be / ~ 3.5 x 10~ 3 Hz, 
with amplitude h — lO -18 ?-^, 1 for = 0.14 and 
h ~ lO -17 ^ 1 for = 0.18. Such signals would be 
easily detectable by the space-based LISA detec- 
tor. 
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Fig. 1 . — Density contours are shown in the x — z 
plane for 4 models with d = 0.2 and N = 3.33. 
The initial maximum density is normalized to 
unity, and the contours are at levels of 0.9, 0.1, 
0.01, and 0.001 

Fig. 2. — 2-D density contours in the equatorial 
plane arc shown for the model with (3 — 0.14 run 
on the cylindrical code. The contour levels are 
0.01 (purple), 0.1 (blue), 0.9 (green), 2 (yellow), 
and 4 (off-white) times the maximum density at 
the initial time, which is normalized to unity. In 
the greyscale version of this figure, the density de- 
creases as the darkness of the shading increases. 
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Table 1: The growth rates, cigenfrcquencies, and 
pattern speeds for the m = 1 and m — 2 modes are 
given in a ring at radius w = 0.32 in the equatorial 
plane. 



Fig. 3. — The growth of the amplitudes \A m \ for 
m = 1 (thick solid line), m = 2 (thin solid line), 
m = 3 (dashed line), and m — 4 (dotted line) is 
shown for the model with initial T/|VF| = 0.14 run 
on the cylindrical hydrocode. These amplitudes 
were calculated in the equatorial plane for a ring 
with radius w = 0.32. 

Fig. 4. — Same as Figure 3, except that the model 
was evolved on the Cartesian hydrocode. 
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